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ABSTRACT 

An harmonic-space maximum-entropy method (MEM) is presented for separating the 
emission from different physical components in all-sky observations by the forthcom- 
ing Planck satellite. The analysis is performed at full Planck resolution, with a pixel 
size of 1.7 arcmin, which corresponds to £ max « 6000. The simulated Planck data 
include emission from the CMB, the kinetic and thermal Sunyaev-Zel'dovich (SZ) ef- 
fects from galaxy clusters, as well as Galactic dust, free-free and synchrotron emission. 
Our simulations also assume homogeneous, uncorrelated pixel noise, although this is 
not a requirement of the method. We find that the MEM technique produces faithful 
reconstructions of the main input components over the whole sky, without the need 
to perform a Galactic cut. The CMB power spectrum is accurately recovered up to 
£ ks 2000. The algorithm is parallelised so that the entire reconstruction can be per- 
formed in ~ 6 hr using 30 R10000 processors on an SGI Origin 2000 supercomputer 
and requires 14 Gb of RAM. 

Key words: methods - data analysis - techniques: image processing - cosmic mi- 
crowave background. 



1 INTRODUCTION 

Two new cosmic microwave background (CMB) satellite 
missions are currently in preparation. The NASA MAP 
satellite is expected to be launched in late 2001, followed 
by the ESA Planck Surveyor in 2007. Both missions will 
provide detailed all-sky maps at a number of observing fre- 
quencies. The main aim of these new satellite projects is to 
obtain an accurate map of CMB anisotropies over the whole 
sky and produce a definitive measurement of the CMB power 
spectrum. This should allow tight constraints to be placed 
on fundamental cosmological parameters and distinguish be- 
tween competing theories of structure formation. 

The maps produced by these satellites will, however, 
contain contributions from various foreground components, 
most notably Galactic dust, free- free and synchrotron emis- 
sion as well as the kinetic and thermal Sunyaev-Zel'dovich 
(SZ) effects from galaxy clusters. In addition, significant con- 
tamination from extragalactic point sources is also expected. 
It is therefore clear that in order to obtain maps of the CMB 
anisotropies alone, it is necessary to separate the emission 
due to these various components. Traditional methods for 
performing the separation include singular-valued decom- 
postion and Wiener filtering (Bouchet, Gispert & Puget 
1996; Tegmark & Efstathiou 1996), although a recent pre- 
liminary application of neural-networks to this problem ap- 
pears promising (Baccigalupi et al. 2000) . A review of tradi- 



tional foreground separation techniques is given by Bouchet 
& Gispert (1999). 

In a previous paper (Hobson et al. 1998); hereafter Pa- 
per I), a separation was performed on simulated Planck Sur- 
veyor observations of a 10 x 10 deg 2 field (see Bouchet et al. 
1997; Gispert & Bouchet 1997), using a non- linear Fourier- 
space maximum-entropy method (MEM) that reduces to 
traditional Wiener filtering in the absence of non-Gaussian 
signals. It was found that faithful reconstructions may be 
produced not only of the CMB anisotropies but also of the 
Galactic components and the thermal SZ effect from mas- 
sive clusters. It was also shown that the MEM technique 
outperformed standard Wiener filtering, particularly in re- 
constructing highly non-Gaussian components such as the 
SZ effects. An application of the MEM component separa- 
tion technique to simulated MAP data of a 10 x 10 deg 2 field 
was presented by Jones, Hobson & Lasenby (1999) 

In Hobson et al. (1999), the basic MEM algorithm was 
extended to identify and remove contamination from extra- 
galactic point sources, and has been further refined by Vielva 
et al. (2001) by combining it with a mexican hat wavelet fil- 
tering technique. 

A great advantage of forthcoming satellite missions, 
however, is the prospect of obtaining all-sky maps at each 
observing frequency, from which one would hope to recon- 
struct maps of the emission from each physical component 
over the whole sky simultaneously, rather than just in small 
patches. Indeed, Prunet et al. (2001) recently presented the 
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first application of the Wiener filter component separation 
algorithm to simulated all-sky data from the MAP experi- 
ment. Their analysis was performed up to ^ max = 512 and in- 
cluded contributions from the three dominant astrophysical 
components expected in the five MAP frequency channels, 
namely CMB and Galactic dust and synchrotron emission. 

It is, however, the Planck mission that provides the 
greatest challenge for performing a component separation. 
Planck will map the whole sky in 10 frequency bands from 
30 to 857 GHz and at angular resolutions ranging from 33 to 
4.5 arcmin. To ensure Nyquist sampling of the sky emission, 
the basic pixel size for Planck reconstructions is ~ 2 arcmin, 
which corresponds to ~ 50 million pixels in each all-sky map 
and maximum multipole ^ max ~ 6000. However, the problem 
is not simply one of computational complexity. Owing to its 
greater frequency coverage and higher angular resolution, 
Planck is sensitive to the full range of astrophysical fore- 
ground outlined above, and thus provides an opportunity 
to map each of these components individually. In particular 
the highly non-Gaussian SZ effects from galaxy clusters are 
of considerable cosmological interest. 

In this paper, we therefore extend the Fourier-space 
MEM algorithm to perform reconstructions on the whole 
celestial sphere. Since the algorithm is inherently computa- 
tionally fast and efficient, and, moreover, can be straight- 
forwardly parallelised to take advantage of existing super- 
computing facilities, we are able to perform reconstructions 
to full Planck resolution. As well as the vastly increased 
computational burden over the analysis presented in Pa- 
per I, when performing all-sky reconstructions we also en- 
counter the problem of highly inhomogeneous emission due 
to Galaxy. The latter leads to foreground maps with a very 
large dynamic range. Nevertheless, we show that presence of 
bright emission from the Galactic plane does not severely af- 
fect the accuracy of the reconstructions obtained, and that it 
is not necessary to impose a Galactic cut on the data prior to 
the component separation analysis. We note, however, that 
a Galactic cut can still be performed in order to prevent con- 
tamination from poorly modelled strong Galactic emission 
regions. 



2 MODEL OF THE MICROWAVE SKY 

To create simulated Planck observations, we must first build 
a plausible model of the emission over the whole sky at each 
Planck frequency. We assume the main contributions to this 
emission are from the primordial CMB, the kinetic and ther- 
mal SZ effects from galaxy clusters, and the Galactic dust, 
synchrotron and free-free components. In this paper, we will 
assume that extragalactic point sources may be removed 
earlier in the analysis using the satellite observations them- 
selves together with existing surveys, or by applying the 
joint MEM and mexican hat wavelet technique discussed 
in Vielva et al. (2001) to the data map at each frequency. 
The generalisation of this joint scheme to all-sky observa- 
tions and using all the frequencies simultaneously will be 
presented in a forthcoming paper. 

Following the notation of Paper I, we choose to work 
in units of equivalent CMB thermodynamic temperature. If 
AI(x) is the fluctuation in the specific intensity at some 
frequency v in the direction x, then the corresponding fluc- 



tuation in the equivalent CMB thermodynamic temperature 
is given by 



AT(x,u) = 



AI(x,v) 



8B(v,T)/dT\ 



T=T 



where B{v, T) is the Planck function and To = 2.726 K is 
the temperature of the CMB (Mather et al. 1994). As in 
paper I, we also assume that the contribution from the pth 
physical component can be factorised into a spatial template 
s p (x) at reference frequency vo and a frequency dependence 
f P (v), so that 



AT(x, v) = At p(*. 

p=i 



(1) 



p=l 



We take as our reference frequency vo = 300 GHz, and nor- 
malise the frequency behaviour so that f p (vo) = 1 for all the 
physical components. 

It is clear that (|l|) represents the somewhat idealised 
case in which the spectral index of each physical compo- 
nent does not vary with direction on the sky. While this 
assumption is valid for the CMB and the two SZ effects, it 
is unlikely to hold for the Galactic components. This was 
not a severe concern in Paper I, since the analysis there was 
restricted to a 10 x 10 deg 2 patch of sky over which the spec- 
tral indices of the Galactic components might reasonably be 
expected not to vary significantly (particularly in regions of 
high Galactic latitude). For an all-sky reconstruction, how- 
ever, the assumption almost certainly does not hold, even 
approximately. Nevertheless, as mentioned in Paper I and 
investigated further in Jones et al. (1999), it is possible for 
the separation algorithm to accommodate varying spectral 
indices by including extra components in the reconstruction. 
For example, if we assume that the frequency dependence 
of the synchrotron emission is of the form / oc , with 
(3 = — 0.7±0.2, we simply include two synchrotron channels, 
one with /3 = —0.5 and one with (5 — —0.9, or even with in- 
termediate values, and afterwards sum over these channels 
to obtain the reconstructed synchrotron map. Alternatively, 
one can consider deviations from the mean spectrum as just 
another template to be recovered with a modified spectral 
behaviour as obtained by linearising the frequency depen- 
dence of the intensity with these deviations (Bouchet et al. 
1996). We will pursue these ideas in a forthcoming paper, 
but for now we consider only the idealised situation repre- 
sented by (Q. 

Each of the six physical components of emission are 
simulated in the HEALPixQ pixelisation scheme (Gorski, 
Hivon & Wandelt 1999) with N sidc = 2048, which corre- 
sponds to ~ 50 million pixels of size 1.7 arcmin, or ^ max = 
3iV s ide — 1 = 6143. Since the highest angular resolution of the 
Planck satellite is 4.5 arcmin, a pixel size of 1.7 arcmin en- 
sures that the fields are Nyquist sampled. The models used 
for each physical component are discussed below. These in- 
put maps are plotted in Fig. [I] at 300 GHz, and the corre- 
sponding power spectra are shown in Fig. H. 



http: / /www. eso.org/kgorski/healpix/ 
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Figure 1. All-sky realisations of the six main physical components contributing to the sky emission at Planck wavelengths. The 
components are primordial CMB, kinetic and thermal SZ effects from clusters and Galactic dust, free-free and synchrotron emission. 
Each map is defined in the HEALPix pixelisation scheme with N sidB = 2048, which corresponds to ~ 50 X 10 6 pixels of size 1.7 arcmin. 
Each map is plotted at 300 GHz in units of fiK. 



2.1 CMB anisotropy 

The primordial CMB anisotropy field is generated over the 
whole sky using the Synfast routine in the HEALPix 
software package, and corresponds a Gaussian realisation 
of a spatially-flat standard inflationary CDM model with 
Q m = 0.35, S7 A = 0.65, Q b h 2 = 0.02 and h = 0.65, for which 
the Ce coefficients were created using Cmbfast (Seljak & 
Zaldarriaga 1996). 



These templates were distributed on the sky according to a 
Poisson distribution and with random orientations. The red- 
shift distribution was chosen to be consistent with the Press- 
Schechter model (Press & Schechter 1974) for a Q m = 0.35, 
Qa = 0.65 cosmology (R. Kneissl, private communication). 
To simulate the kinetic SZ effect, the cluster radial velocities 
are assumed to be Gaussian distributed with a dispersion of 
400 km s" 1 at z = 0. 



2.2 Thermal and kinetic SZ effects 

The thermal SZ effects from individual clusters were simu- 
lated using the gas dynamics code of Eke, Navarro & Frenk 
(1998), from which 210 cluster templates were obtained. 



2.3 Galactic dust emission 

The Galactic dust contribution is modelled using the 
DIRBE-1RAS 100-^m dust map (Schlegel, Finkbeiner & 
Davis 1998). The angular resolution of the map is around 5 
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arcmin, which is comparable with Planck resolution in the 
highest frequency channels. The emission at each Planck fre- 
quency can be predicted by extrapolating the 100-/jm flux 
assuming a one-component dust model with a temperature 
Tdust = 18 K and dust emissivity = 2. The colour correc- 
tion factor for the DIRBE 100-/im filter was also taken into 
account; (see Finkbeiner, Davis & Schlegel 1999). 



2.4 Galactic synchrotron emission 

The basic template of the synchrotron emission is the de- 
striped version of the 408 MHz Haslam survey (Haslam et al. 
1982), to which additional structure was artificially added 
at sub-degree angular scales by extrapolating the angular 
power spectrum of the Haslam survey as Ci oc The 
synchrotron map at 300 GHz is obtained by extrapolating 
the 408 MHz survey using an all-sky map of the spectral in- 
dex constructed by combining the low frequency surveys at 
408 MHz, 1420 MHz (Reich & Reich 1986) and 2326 MHz 
(Jonas, Baart & Nicholson 1998) and padding the unob- 
served area around the South pole with the mean spectral 
index, at a resolution (FWHM) of 10 degrees (G. Giardino, 
private communication). To predict the synchrotron contri- 
bution at the Planck observing frequencies, a constant in- 
tensity spectral index (5 — —0.9 was assumed. 



2.5 Galactic free- free emission 

Reliable maps of Galactic free- free emission are not currently 
available, although experiments such as the Ha Sky Sur- 
ve yQ (Gaustad et al. 2001) and the Wisconsin Ha Mapper 
- WHAM | (Haffner 2001) should soon provide maps that 
could be used as templates. The resolution of the surveys is 
about 0.8 arcmin for Ha Sky Survey and about 1 degree for 
WHAM. For the time being, however, we create a free-free 
template based on the DIRBE/IRAS dust map, of which 60 
per cent of the emission is a dust-correlated component and 
40 per cent is uncorrelated. The uncorrelated component is 
simply a dust map flipped north-south in Galactic coordi- 
nates. The spectral index of the free-free emission is assumed 
to be P = —0.16 and the normalisation is that suggested by 
Bouchet & Gispert (1999). 



3 SIMULATED PLANCK OBSERVATIONS 

If we observe the sky at n/ observing frequencies then, in 
any given direction x, we obtain a data vector of length 
71/ that contains the observed temperature fluctuation in 
this direction at each observing frequency plus instrumental 
noise. In order to relate this data vector to the emission from 
each physical component it is useful to introduce the rif xn c 
frequency response matrix with components defined by 

F vp = / tu(v')fp(y')&v , 
Jo 



2 http:/ /www. swarthmorc.edu/Homc/Ncws/Astronomy/ 

3 http:/ /www. astro. wisc.edu/wham/ 
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Figure 2. Power spectra of the input maps shown in Fig. hi 



where t v {v') is the frequency response (or transmission) of 
the vth. frequency channel. Assuming that the satellite ob- 
serving beam in each channel is spatially- invariant, we may 
write the beam-smoothing as a convolution and the ^th com- 
ponent of the data vector in the direction x is then given 
by 

B v (x ■ x) F vp s p (x ) dCl' + e v {x) (2) 
-*~ P =i 

where B v is the beam profile for the z/th frequency channel. 
The beam profile at each frequency is assumed to be Gaus- 
sian. The e v (x) term represents the instrumental noise on 
the observations in the vth. channel, which for simplicity we 
assume to be uncorrelated Gaussian noise with a fixed rms 
over the whole sky. The assumed observational parameters of 
the Planck satellite are given in Table 0. Fig. | shows the rms 
equivalent thermodynamic temperature fluctuation at each 
Planck frequency due to each physical component, in the sky 
region lying between Galactic latitudes 65° < b < 75°. The 
rms instrumental noise per pixel in each frequency channel 
is also plotted. 

If the beam profile at each frequency is spatially- 
invariant and circularly-symmetric, it is convenient to work 
in terms of spherical harmonic coefficients. We now consider 
the statistical properties of these coefficients. 

3.1 All-sky observations 

If each data map d v (x) is defined over the whole sky, then on 
performing a spherical harmonic transform the convolution 
in (|^) becomes a simple multiplication. Thus, we obtain 

Tic 

p=l 

where we have adopted the usual notation for spherical 
harmonic coefficients ft m — dfl Y e * m (x)f(x) in which 

Y( m (x) is a standard spherical harmonic function. Thus, 
are the spherical harmonic coefficients of the ^th frequency 
map. These are related to the spherical harmonic coefficients 
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Table 1. The assumed observational parameters for the Planck Surveyor satellite. Angular resolution is quoted as FWHM for a Gaussian 
beam. Sensitivities are quoted per beam FWHM for 12 months of observation. 



Low Frequency Instrument High Frequency Instrument 



Central frequency (GHz): 


30 


44 


70 


100 


100 


143 


217 


353 


545 


857 


Fractional bandwidth (Au/u): 


0.2 


0.2 


0.2 


0.2 


0.37 


0.37 


0.37 


0.37 


0.37 


0.37 


Transmission: 


1.0 


1.0 


1.0 


1.0 


0.3 


0.3 


0.3 


0.3 


0.3 


0.3 


Angular resolution (arcmin): 


.33 


23 


14 


10 


10.6 


7.1 


4.9 


4.5 


4.5 


4.5 


AT sensitivity (//K): 


4.4 


6.5 


9.8 


11.7 


4.9 


5.7 


12.5 


40.9 


392 


12621 
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Figure 3. The rms equivalent thermodynamic temperature fluc- 
tuation per pixel at each Planck observing frequency in the region 
of sky with Galactic latitude 65° < b < 75°. The rms instrumen- 
tal noise per pixel in each frequency channel is also plotted. The 
pixel size is 1.7 arcmin. 



a lm °f the pth physical component via the response matrix 
for the observations R^" = F vp , where are the 

harmonic coefficients of the ^th observing beam. Finally, 
e trn i s * n e instrumental noise on the harmonic mode (£,m). 
It is important to note that (^) is satisfied at each spher- 
ical harmonic mode (£, m) independently. Thus, in matrix 
notation, at each mode we have 

dim = Re^tm + Sim, (4) 

where de m , ae m and ei m are column vectors containing n/, 
n c and nj complex components respectively, and the re- 
sponse matrix has dimensions n/ x n c . 

Although the data at a given (£, m) mode do not de- 
pend on the harmonics of the physical components at other 
modes, the a priori covariance structures of the harmonic co- 
efficients of the components azl and the instrumental noise 
Cpfl do correlate different modes. The sky emission (from the 
Galactic components) is anisotropic owing to the presence 
of pronounced emission from the Galactic plane. Further- 
more, scanning strategies with non-uniform coverage of the 
sky lead to spatial variations in the noise rms per pixel, and 
so the instrumental noise field is also anisotropic. This re- 
sults in a priori correlations between different (t,m) modes, 
so that in general we have 

(a^ m a^, m ,) = Q m ,£' m /, (5) 



{€l m e\, m ,) — ^tm,l'm', (6) 

where we have defined the set of ensemble-average signal 
and noise covariance matrices C^ m ,f m ' and Ne m ,i'm'- Each 
matrix Qtm,i'm' has dimensions n c x n c and contains the a 
priori covariances of the harmonic coefficients of the different 
physical components at the modes (£,m) and (£' ,m'). Each 
matrix N fm f / m / has dimensions nfXnf and contains the co- 
variances between the modes of the instrumental noise har- 
monics at the different observing frequencies. It should also 
be remembered that the emission from components other 
than the CMB is highly non-Gaussian and so these second- 
order statistics will not, in general, provide a full statistical 
description of the sky emission, particularly in the Galactic 
plane. In principle, we can also include any prior informa- 
tion about noise correlation properties in the noise covari- 
ance matrix N. A typical example is non-homogeneous pixel 
noise in each frequency channel. 

3.2 Galactic cuts 

Since Planck is primarily a CMB instrument, 'component 
separation' is often taken to mean 'foreground removal'. In 
particular, mapping of Galactic emission is regarded only as 
a secondary goal of the mission. When performing an all- 
sky foreground removal, one might therefore consider first 
removing the pronounced emission from the Galactic plane 
by imposing a Galactic cut on each individual frequency 
map. The new input data maps are then given by 

d v (x) = W(x)d v (x), (7) 

where the window W(x) is zero in the Galactic cut and 
equals unity outside it (although more general windows are 
possible). As indicated in (^), the window W is usually cho- 
sen to be the same for each observing frequency, and is often 
taken to be symmetric about the Galactic equator (for ease 
of computation, as discussed below). Indeed, this approach 
was adopted by Prunet et al. (2001). 

The imposition of a Galactic cut does, however, produce 
some troublesome effects. Firstly, the spherical harmonic ba- 
sis Yi m (x) is no longer orthonormal over the cut-sphere. For 
the Wiener filter approach, in particular, this leads to in- 
consistencies between performing the component separation 
in real space and harmonic space (see Prunet et al.). Sec- 
ondly, the Galactic cut leads to additional correlations be- 
tween harmonic coefficients at different (£, m) modes. For a 
function f(x) defined on the sphere, the harmonic modes 
after the Galactic cut / = Wf are given by 

hm = / dttYi m (x)W(x)f(x) = W im,i'm'fe'm', 
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where the coupling matrix is 



W, 



lm,,t' ' m' 



dnw(x)Y; m (x)Y e , m ,(i 



The problem of loss of orthonormality is usually solved 
by performing the entire analysis in terms of a new basis 
Y' tra {x) that is orthonormal on the cut-sky, such that 



dQ W (x)Y'i m (x)Yl, m , (x) = 5 U 'S V 



For ^max < 50, this basis can be constructed by performing 
a Cholesky decomposition of the coupling matrix (Gorski 
1994). For higher £ ma x, a singular value decomposition 
(SVD) is required (Tegmark 1997; Mortlock, Challinor & 
Hobson 2001). We note that, for simulated MAP data with 
^max = 512, Prunet et al. find the computational and mem- 
ory costs of computing the new orthonormal basis to be rea- 
sonable, but the corresponding computing requirements for 
Planck will be somewhat larger. Perhaps more importantly, 
the computation of the new orthonormal basis is entirely 
unfeasible for either MAP or Planck unless the edges of the 
Galactic cut are chosen to be lines of constant latitude. It 
should also be remembered that for the new basis functions 
the subscripts I and m no longer have their usual physical 
meanings. Nevertheless, for large £ ma x (such as for MAP or 
Planck) the new basis functions Y' tm (x) are very close to 
zero in the cut, and so the coefficients of these modes have 
the desirable property of being insensitive to emission from 
the cut region. 

Unfortunately, the construction of the new orthonor- 
mal basis does not solve the problem of mode coupling. The 
new basis can be expressed as a linear combination of the 
standard spherical harmonic basis functions, 



Yf„ 



(#) = £ Otm,t'm'Yl'm'(&)> 



t' ,m' 



where matrix elements 0< m /' m ' are given by the 'overlap' 
integrals of the new and old basis functions, 



O, 



m.t' m' 



AO. Y'\ m (x)Y llm ,{x). 



If we expand the sky emission from the physical components 
in terms of our new basis, we thus obtain 



im.l' m 1 a l' m' > 



(8) 



and a similar expression holds for the coefficients e' lm of the 
instrumental noise field at each observing frequency. Hence, 
the a priori correlation structure in the new basis is given 
by 



J (m d I'm 1 



— £'in 

= N'„ 



(9) 
(10) 



where explicit expressions for e , m , and H' lr 
ily obtained from (^), (|B|) and (jq). In general, coupling still 
persists betweem different (£,m) modes. 

An additional problematic feature of performing a 
Galactic cut is that the operation of beam convolution is 
no longer diagonal in the new orthonormal basis. From (^|) 
and (H), we have 



d' 



p=l V ' jn' 



(up) 



» 



/(") 



(11) 



which, in general, cannot be written in the form 

+ e lm . (12) 

and this leads to further coupling between different (I, m) 
modes. Thus, in general, the operation of beam convolution 
is not merely non-diagonal, but is not strictly defined in a 
formal sense. We note, however, that we can write ( [tl| ) in 
this form in the special case where the response matrix R is 
not a function of £. This corresponds to no beam convolu- 
tion. 



3.3 Incomplete sky coverage 

Similar considerations to the above also apply to incom- 
plete sky coverage resulting from, for example, the scanning 
strategy of the satellite or instrument failure. In this case, 
however, W can be highly irregular in shape and different 
for each observing frequency. If the affected regions are not 
too large one could define a single window W that encom- 
passes the unobserved regions at all frequencies, but this 
would probably not be viable in practice. Moreover, it is 
unlikely that one could find an adequate Galactic cut with 
edges at constant latitude, and so it would be computation- 
ally impossible to construct an new orthonormal basis on 
the cut-sky. 



3.4 Choice of basis functions 

Following the above discussion, we have the choice of per- 
forming the calculation in terms of the standard spherical 
harmonic basis functions Yi m or some other basis Y[ m that 
is orthonormal on the part of the sky outside some con- 
stant latitude Galactic cut. In both bases, there is coupling 
of modes through their a priori covariance structure, as de- 
fined by (^-^) or (p|-|ici|) respectively. Thus, in either case, a 
general analysis of the data cannot be performed 'mode-by- 
mode', but requires one to estimate the full 'signal' vector 
a (or a') simultaneously from the full data vector d (or d'); 
these vectors can be regarded simply as the concatenation of 
the vectors a^ m and Ag m respectively for all possible values 
of I and m (or their primed counterparts) . 

It is also worth noting, however, that in the cut-sky 
basis one must assume that (|l2|) holds, which is not valid 
unless there is no beam convolution. This assumption is, of 
course, unnecessary in the standard harmonic basis, since 



(|) holds exactly for all-sky observations. The assumption 
(121) will, in general, lead to reconstruction artefacts which 
are most pronounced near the edges of the Galactic cut. 

Owing to the above difficulties of working in the cut-sky 
basis, and given that we are interested in recovering maps of 
the Galactic components as well as the CMB, in this paper 
we choose not to impose a Galactic cut, and instead attempt 
to reconstruct the physical components over the whole sky. 
Thus we perform our calculations in the standard spherical 
harmonic basis. 
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4 COMPONENT SEPARATION 

The aim of any component separation algorithm is to use 
the data maps d v {x) at each frequency to obtain an esti- 
mated maps s p (x) of the emission from each physical com- 
ponent. Typical methods include singular-valued decompo- 
sition, Wiener filtering or the maximum-entropy method, all 
of which can be viewed in Bayesian context (see Paper I). 

4.1 Harmonic-space maximum-entropy method 

A complete description of the 'mode-by-mode' Fourier-space 
MEM component separation algorithm as applied to small 
patches of sky is given in Paper I. The main difference in 
applying the technique to all-sky reconstructions is that 
Fourier transforms must be replaced by spherical harmonic 
transforms. Aside from this modification, the algorithm re- 
mains basically unchanged, and so we give only a brief out- 
line here. 

Let us consider the standard (unprimed) spherical har- 
monic basis, although the following discussion is equally 
valid for the primed basis described above. Using Bayes' the- 
orem, the estimator a of the signal vector is usually taken 
to be that which maximises the posterior probability 

Pr(a|d) oc Pr(d|a) Pr(a) 

where Pr(d|a) is the likelihood of obtaining the data given 
a particular signal vector and Pr(a) is the prior probability 
that codifies our expectations about the signal vector be- 
fore acquiring any data. Indeed, it is the specification of the 
prior alone that differentiates between different component 
separation algorithms such as singular-value decomposition, 
Wiener filtering and the maximum-entropy method (see Pa- 
per I). 

Unfortunately, the lengths of the (complex) data and 
signal vectors are ~ n./^max and ~ nc'maj, and so for MAP 
or Planck observations the dimensionality of the numerical 
maximisation problem renders it computationally unfeasi- 
ble. It is therefore necessary to make some approximation. 
The most straightforward approximation is to neglect any 
coupling between different (I, m) modes, and perform the re- 
construction independently 'mode-by-mode'. Although this 
may appear an extreme approximation, we will see in Sec- 
tion |E| that it still produces excellent reconstructions. Taking 
the modes to be independent corresponds to assuming that 
the likelihood and prior probability distributions factorise, 
such that 

Pr(d|a) = ][[Pr(d, m |af m ), (13) 

t,m 

Pr(a) = n Pr ( a ^- ( 14 ) 

This offers an enormous computational advantage, since one 
can maximise the posterior probability 

Pr(a £m |df m ) oc Pr(d 

£m|^m) Pr(a^ m ) 

at each mode separately. Thus one replaces a single numer- 
ical maximisation over a parameter space with 2n c Z^ lax real 
dimensions by Z^ ax separate minimsations with 2n c real di- 
mensions. 

The factorisations (^) and ( |l^ ) imply an assumed a 



priori covariance structure of the form 
(af m a|, m ,) = Ci m Sei'S mm i , 

Although not required by the mathematical formalism, we 
have made the additional simplifying assumption that Ct m 
and Nt m do not depend on m, so that 

(a^ TO a f , TO ,) = CeS ee iS mrn i , (15) 

(e« m e], TO ,) = NtSu'Smm'- (16) 

In the standard spherical harmonic basis, this approxi- 
mation is equivalent to assuming that the emission from 
each physical component and the instrumental noise are all 
isotropic random fields on the celestial sphere. The matri- 
ces Ce and contain the (cross) power spectra of, respec- 
tively, the physical components at vq and the instrumental 
noise at each observing frequency. Unfortunately, the analo- 
gous approximation in the cut-sky basis does not have such 
a straightforward interpretation. 

At each (£, m) mode, any a priori correlation informa- 
tion is contained in the matrix Ct. As explained in Paper I, 
this information is most easily incorporated into the recon- 
struction algorithm by first defining a 'hidden' vector he m 
at each mode, that is related to the vector a^ m by 

ai m = Llhlm- (17) 

The n c x n c lower triangular matrix is obtained by per- 
forming the Cholesky decomposition Ce = L^Lj. The com- 
ponents of the hidden vector thus have the useful property 
of being a priori uncorrelated with unit variance, since then 

(ae m a\ m ) = (Uh^ ro h] m L]) = U{h tm h\j\-\ = C e . 

The reconstruction is then performed entirely in terms of 
hem and the corresponding reconstructed signal vector is 
subsequently found using (|l^). 

Assuming the instrumental noise to be Gaussian, the 
likelihood function at each mode is given by 

Pr(d oc exp [-x 2 (hf 

m ) \ , 

where \ 2 is the standard misfit statistic given by 
X 2 (hem) = (d^ m — RzLebe m ) IM^ 1 (d^ m — R^L^h^ m ) 
For the prior Pr(h^ m ), we assume the entropic form 
Pr(h em ) oc exp[aS(h em , m)] 

where S(he m , m) is the cross entropy of the complex vectors 
hem and m, where m is a model vector to which hi m defaults 
in absence of data. The form of the cross entropy for complex 
images is given in Paper I. Thus, maximising the posterior 
probability Pr(h^ m |d^ m ) with respect to he m is equivalent to 
minimising the function 

$(hft») = X 2 Om) - aS{h lm , m). (18) 

As explained in Paper I, the regularising parameter a 
can be determined in a Bayesian manner by treating it sim- 
ply as another parameter in the hypothesis space. Indeed, 
by assuming a Gaussian form for the posterior probability 
distribution near its maximum, one can derive a closed-form 
expression for the optimal value of a (see Paper I). Since the 
calculation is performed at each (£, m) mode independently, 
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one could in principle obtain the optimal value of a sepa- 
rately for each mode. However, since the dimensionality of 
the parameter space at each mode is only 2n c -dimensional 
(i.e. 12-dimensional in our case), the Gaussian approxima- 
tion to the posterior distribution is occasionally too inac- 
curate to determine a properly. Nevertheless, we would not 
expect the optimal level of regularisation to depend strongly 
on m, but to vary considerably with I. Thus we assume the 
optimal value of a to be the same for all modes at the same 
t. Taking these modes together, the dimensionality of the 
parameter space it then sufficiently large that the Gaussian 
approximation to the posterior distribution is reasonably ac- 
curate, and the optimal value of a is well-determined. Thus, 
the harmonic-space MEM algorithm employs optimal, scale- 
dependent regularisation, which allows for a more accurate 
reconstruction. 

The numerical minimisation of ( |l8| ) can be performed 
using a variable metric minimiser (Press et al. 1994). As a 
result of the independence of each (£,m) mode, paralleli- 
sation of the separation algorithm is straightforward. Full- 
resolution Planck simulations in HEALPix with iV s ide = 
2048 have a maximum multipole ^ max = 3iV s id e — 1 = 6143. 
The entire component separation can be performed in ~ 6 
hr using 30 R10000 processors on the ' COSMOS ' SGI Origin 
2000 supercomputer and required only 14 Gb of RAM. 

As shown in Paper I, one can also estimate the covari- 
ance matrix of the reconstruction errors on the harmonic 
modes at m - Making a Gaussian approximation to the pos- 
terior distribution at each mode, we find 

((aim — aim) (aim — aim)'") = L^H^Lj, (19) 

where H^ m = VV$(hf,„) is the Hessian matrix of the poste- 
rior distribution at its peak hf m . For full-resolution Planck 
data, the calculation of the errors on the component separa- 
tion required just 15 mins additional computation time on 
Cosmos. 



5 APPLICATION TO SIMULATED DATA 

We now apply the separation algorithm outlined above to 
the simulated Planck Surveyor observations described in 
Section Owing to the limited resolution of our foreground 
maps, we only consider multipoles up to ^ ma x = 4096, which 
is the Nyquist sampling limit for pixels of size 1.7 arcmin. 
This reduces both the computation time and memory re- 
quired by about a factor of two, as compared with the case 
U = 6143. 

Clearly, the quality of the reconstructions will depend 
on our prior knowledge of the input components. Following 
Paper I, we assume that the frequency spectrum behaviour 
of the components is accurately known. Nevertheless, as 
mentioned in the introduction, the relaxation of this assump- 
tion has investigated by Jones et al. (1999), who showed 
that it is possible to accommodate mild uncertainties in fre- 
quency behaviour and spatial variation of spectral indices. 
Prior power information may be incorporated into the algo- 
rithm through the signal covariance matrix Ci given in (|l5|). 
Similarly, any knowledge of the correlation structure of the 
instrumental noise can be included via the noise covariance 
matrix Nf. In our simulated data, however, we assume homo- 



geneous, uncorrelated instrumental noise in each frequency 
channel. 



5.1 Partial prior covariance information 

As an illustration, we begin by assuming knowledge of the 
azimuthally-averaged power spectra of each input compo- 
nents, as shown in Fig. ^, except that of the CMB. Since the 
accurate determination of the positions and heights of the 
acoustic peaks in the CMB spectrum is of such fundamental 
importance, we choose not to prejudice our reconstruction 
by supplying the component separation algorithm with the 
power spectrum of the input CMB map. Thus, as our as- 
sumed CMB power spectrum, we take a simple functional 
form which matches the overall shape of the true spectrum, 
but contains no acoustic peaks (this assumed power spec- 
trum is plotted as the dotted line in the first panel of Fig.^|). 

The inclusion of the azimuthally-averaged input power 
spectra for the remaining components is reasonable, and 
would be possible in practice. As discussed in section H, the 
maps of the kinetic and thermal SZ effect have been pro- 
duced with some care, and we would expect them to pro- 
vide a reasonably accurate representation of the true power 
spectra for these two components. Alternatively, one could 
use an average power spectrum for each SZ component ob- 
tained from a large number of simulated realisations, which 
would be statistically more accurate. For the Galactic com- 
ponents, the dust and synchrotron maps are simply taken 
from existing observations (albeit with some added structure 
on small scales) . Unfortunately, considerable uncertainty re- 
mains for Galactic emission in the free-free. Nevertheless, 
free-free observations should be available before the Planck 
data are analysed, so that an approximate free-free power 
spectrum could be calculated. As discussed in section bl for 
the time being, we create a free-free template based on the 
DIRBE/IRAS dust map, of which 60 per cent of the emission 
is a dust-correlated component and 40 per cent is uncorre- 
lated, with the uncorrelated component being simply the 
dust map flipped north-south in Galactic coordinates. 

We note, however, that, since we are assuming knowl- 
edge of only the power spectra of each component, the com- 
ponent separation algorithm has no a priori information con- 
cerning cross correlations between any of the components. In 
practice, a significant amount of cross correlation informa- 
tion should be available to the Planck mission from earlier 
observations, and so our assumed level of prior knowledge 
is, in this respect, mildly pessimistic. 



5.1.1 Reconstructed and residuals maps 

The component separation algorithm produces an estimate 
a£ m of the harmonic coefficients of the physical components 
at each mode. Let be the true harmonic coefficient of the 
pth physical component at the mode (£,m). This coefficient 
can be written 

(p) _ a (p) , g (p) (on) 

where a?l is the estimate of the mode obtained from the 
component separation algorithm and 8a^ is the residual 
on this mode. The reconstructed map of the pth physical 
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component is given simply by performing an inverse spheri- 
cal harmonic transform on the reconstructed coefficients a!f^ 
using the HEALPix package, ft is also of interest to con- 
struct a map of the reconstruction residuals for each physical 
component. This can be calculated in two equivalent ways. 
One can either simply subtract the reconstructed map from 
the true input map, or one can perform an inverse spherical 
harmonic transform of the residuals 5af^. We have verified 
that these two methods lead to residuals maps that are con- 
sistent to within the machine precision. 

In Fig^, we plot the reconstructions and residuals of the 
extragalactic components, namely the CMB and the kinetic 
and thermal SZ effects. We see that the CMB reconstruc- 
tion is a faithful respresentation of the input map shown in 
Fig. jjj. This is confirmed by the CMB residuals map, which 
contains very little structure. Indeed, only a thin band of 
residual emission remains in the CMB reconstruction in the 
very centre of the Galactic plane. Thus, the CMB map has 
been faithfully recovered over the whole sky, without the 
need to perform a Galactic cut prior to the component sep- 
aration. It is also worth noting that the true input CMB 
map was not smoothed before calculating the reconstruc- 
tion residuals. FigJ^ shows the distribution of pixel values in 
the residuals map. This distribution is well approximated by 
a Gaussian centred on zero with an rms of 12.3 fiK, except 
for low-level wings corresponding to the very few pixels in 
the Galactic plane that contain residual Galactic contami- 
nation. 

In the absence of a priori cross-correlation information, 
we see that the recovery of the kinetic SZ effect is very poor, 
confirming the findings of Paper I. This result is not surpris- 
ing, since the kinetic SZ effect has the same frequency spec- 
trum as the primordial CMB but, as shown in Figs ^| and 
|§| its rms signal is 3 orders of magnitude smaller than the 
CMB component. However, the thermal SZ effect has been 
successfully recovered, particularly in the brighter clusters. 
The numbers of clusters recovered and accuracy achieved, 
together with the resulting cosmological information that 
could be obtained, will be discussed in detail in a forthcom- 
ing paper. 

In Fig.^| we plot the reconstructions and residuals of 
the Galactic components, namely dust, free-free and syn- 
chrotron emission. Once again we see that each component 
has been successfully recovered over an extremely wide dy- 
namic range in pixel values. In particular, the presence of 
strong emission from the Galactic plane clearly invalidates 
the assumption of uncorrelated harmonic modes, as dis- 
cussed in section ^. Nevertheless, this approximation does 
not appear to have had a dramatic effect on the quality of 
the reconstructions. The accuracy of the reconstructions is 
confirmed by the corresponding residuals maps, which are 
mostly featureless apart from some slight residual emission 
in the central Galactic plane. The histograms of the pixel 
values in these residual maps are shown in Fig. ^. As for the 
CMB, each histogram is well approximated by a Gaussian 
centred on zero, expect for low-level wings corresponding to 
the few pixels in the Galactic plane that contain residual 
contamination from other components. The parameters of 
the best-fit Gaussian to each of the residuals histograms in 
Figs P & § are given in Table [| We note that the low val- 
ues of a for the kinetic SZ effect and Galactic synchrotron 
emission do not indicate that they are the most accurately 



Table 2. The mean fi and standard deviation a (in fiK) of the 
best-fit Gaussian to each of the residuals histograms shown in 
Figs|&§. 



Component 


/' 


<7 


CMB 


0.11 


12.30 


Kinetic SZ 


0.02 


0.12 


Thermal SZ 


0.62 


4.30 


Dust 


0.14 


16.70 


Free- Free 


0.08 


2.40 


Synchrotron 


3 x 10" 4 


0.11 



reconstructed, but simply that the level of emission in these 
two components is very low at 300 GHz. 



5.1.2 Power spectra of reconstructed and residuals maps 

In the left-hand columns of Figs |^ & ^ , we plot the power 
spectrum of the reconstructed map for each physical com- 
ponent (dotted line), together with the power spectrum of 
the true input map (solid line). For the pth physical compo- 
nent, the power spectrum of the reconstructed map is given 
simply by 



C 



2£+l 

m=-£ 



(21) 



In the right-hand column of Figs Q & 0, we plot the power 
spectrum of the residuals map for each component (solid 
line), which is given by 



sc: 



21- 



E 



\Sa 



(p)|2 



(22) 



Of course, in practice the true residuals Sai m would not be 
available. Nevertheless, we can estimate the residuals power 
spectrum by calculating the Hessian matrix V\i m of the pos- 
terior distribution at its peak. As discussed in section ^, 
this provides an estimate of the standard error on the re- 
constructed harmonic modes. From (|l9|), we see that for the 
pth physical component |Aa^j 2 = (|^a^| 2 ) is given by the 
pth diagonal entry of the matrix L^H^L]. Thus, our esti- 
mator of the residuals power spectrum is simply 



SC. 



(p) 



2f- 



E 



»|2 



(23) 



; (p) 



which is unbiassed in the sense that SCt' = (SC { e p) ), where 
the averages are over realisations of the instrumental noise 
and maps of the physical components. The estimated resid- 
uals power spectrum for each component is plotted as the 
dotted line in the right-hand columns of Figs || & ^. 

Fig. |^ shows the power spectra for the extragalactic 
components, namely the CMB and kinetic and thermal SZ 
effects. For the CMB, in addition to the power spectra of 
the input and reconstructed maps, we also plot the sim- 
ple a priori form of the power spectrum assumed by the 
component separation algorithm. This assumed power spec- 
trum contains no acoustic peak structure, so as not to prej- 
udice the reconstruction. We see, however, that the power 
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CMB restored CHE residuals 




-i.ee ^^^^^^^^^^^h -^^^^^^^m +25. ee -12.se ^^^^^^^^^^^m ^^^^^^^h +12. se 

Figure 4. Reconstructions of the extragalactic components, namely the primordial CMB, the kinetic SZ effect and the thermal SZ effect 
(left-hand column). Also plotted are maps of the corresponding reconstruction residuals for each component (right-hand column). All 
maps are plotted in units of fiK. 
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Figure 5. Histograms of the residuals maps for primordial CMB (left), the kinetic SZ effect (centre) and the thermal SZ effect (right). 
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Figure 6. Reconstructions of the Galactic components, namely the dust, free-free and synchrotron emission (left-hand column). Also 
plotted are maps of the corresponding reconstruction residuals for each component (right-hand column). 




Figure 7. Histograms of the residuals maps for emission from dust (left), free- free (centre) and synchrotron (right). 
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CMB reconstruction CMB residuals 
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Figure 8. Power spectra of the reconstructed maps (dotted line) and input maps (solid line) of the extragalactic components, namely 
CMB, kinetic SZ and thermal SZ (left-hand column). The smooth dotted line in the CMB panel is the assumed input CMB power 
spectrum supplied to the algorithm. Also plotted are the power spectra of the residuals maps (solid line) and the predicted residuals 
power spectra (dotted line) for each component (right-hand column). 



spectrum of the reconstructed map matches that of the in- 
put map up to I ~ 2000. In particular, the first 5 acoustic 
peaks have been accurately recovered, and there is some ev- 
idence of the presence of a 6th peak in the recovered power 
spectrum, although the amplitude is significantly underes- 
timated. Indeed, beyond I w 2000, the power spectrum of 
the reconstructed map consistently underestimates that of 
the true map. This is effect is due to the fact that (J2TJ) is a 
biassed estimator of the true power spectrum. Nevetheless, 
as discussed below, we can easily construct an unbiassed es- 
timator. From Fig. H, however, we see that our estimator for 



the power spectrum of the CMB residuals is already unbi- 
assed, as we expected. Thus the MEM algorithm provides a 
good estimate of the statistics of the reconstruction errors. 

The power spectra of the reconstructed kinetic SZ map 
severely underestimates the true power spectrum, as we 
might have expected. Indeed, one can see that the 'recov- 
ered' kinetic SZ map has a power spectrum which clearly 
mirrors that of the CMB, which suggests that any structure 
present in the kinetic SZ 'reconstruction' is dominated by 
contamination by primordial CMB. Neverthelss, once again 
the estimated residuals power spectrum is a good approxi- 



All-sky component separation for the Planck mission 13 



Dust ncconslvuclion 




I a i <n man 



Mullpde number, r 



Free— fLee reconstruction. 




i 10 ■« 



Mullpde nurrber. i 
Figure 9. As in Fig. H, but for the Galactic 



mation to the power spectrum of the true residuals map. In 
practice this would provide a robust indication of the poor 
quality of the kinetic SZ reconstruction. 

For the thermal SZ effect, the power spectrum of the 
reconstructed map matches that of the input map up to I w 
1000, beyond which the true spectrum is underestimated. 
This is again is result of (|l|) being an biassed estimator. 
We see, however, that the predicted and actual residuals 
power spectra are in close agreement for all I values. 

The power spectra of the reconstructed maps and resid- 
uals of the Galactic components are shown in Fig. ^[ Quali- 
tatively, we see similar behaviour to the observed for the ex- 
tragalactic components, with the exception that the power 
spectrum of the input dust map is accurately recovered for 
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components dust, free-free and synchrotron. 



all I. For the free- free and synchrotron components, the un- 
derestimation of the true power spectrum is again due to 
the fact that is a biassed estimator. In each case, we 
again note that the estimated and actual residuals power 
spectrum agree well. The accuracy of this error estimation 
enables one to define an unbiassed estimator of the true in- 
put power spectrum. 



5.1.3 An unbiased power spectrum estimator 

We see from the left-hand columns of Figs ^ & ^, that the 
power spectra of the reconstructed maps generally underes- 
timate the true input power spectra at large I. Indeed, it is 
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easy to show that the estimator is biassed. From (p(i|), 



we have af^ 



» 



5a^, and on inserting this expression 



■ I t, II 6 £. //(. 1,1 lb' — *■ 

into (Elf), we see that the expectation value of this estimator 
is given by 



(p)\ 



where the ensemble average is again taken over realisations 
of the instrumental noise and of the physical components. 
Using (^) again, we can write 



,(P)|2\ 



(24) 



If we assume that the reconstruction errors on each mode are 
not correlated with the reconstructed values (which is rea- 
sonable given the featureless nature of the residuals maps), 
the first term on the right-hand side of (j24|) disappears, and 
so we obtain 



2£+l 

m=-t 

(Cf >) - (5C™), 

{ t } V «&) ) 



(|a^| 2 )-(l^| 2 >, 



(25) 



Thus, ( |2l[ ) is a biassed estimator, which underestimates the 
true power spectrum of the pth component on angular scales 
where the residuals are large. 

By analogy with the approach developed by Bouchet & 
Gispert (1999) and Knox (1995), we denote the quantity in 
brackets in ( p5j ) as the quality factor . We see immedi- 
ately that an unbiassed estimator of the power spectrum is 
given by simply by 



■ — (p) 

5CT 



(26) 
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Figure 10. The power spectrum of the input CMB map (dashed 
line) compared with the unbiassed estimator (Eq) of the power 
spectrum obtained from the CMB reconstructionfgrey solid line). 
Also plotted are the 1-tr error limits. The smooth dotted line 
shows the power spectrum of the CMB assumed as our prior in 
the component separation process. 



Table 3. The mean fi and standard deviation a (in fiK) of the 
best-fit Gaussian to the residuals histograms for full a priori cor- 
relation information. 



Component 


H 


(j 


CMB 


0.09 


12.18 


Kinetic SZ 


0.02 


0.09 


Thermal SZ 


0.60 


4.50 


Dust 


0.14 


16.70 


Free- Free 


0.05 


1.58 


Synchrotron 


5 x 10" 4 


0.10 



—~(p) I — 
where SCe is given by (G2 

is given by 



V 



2£+l 



The variance of this estimator 



(27) 



which illustrates the uncertainty added by the foreground 
removal to cosmic variance. 

As an illustration, in Fig. ^ we plot the unbiassed es- 
timator of the CMB power spectrum against the true power 
spectrum of the input map. We see from the plot that the 
estimated power spectrum is now fully consistent with the 
input power spectrum out to large £ values. In particular, 
the heights and positions of the 5th and 6th acoustic peaks 
are now accurately recovered. Beyond the 6th peak, how- 
ever, there is no evidence for further oscillatory structure in 
the power spectrum. At high £, the overall level of power 
is recovered correctly, but this is simply a result of our as- 
sumed power spectrum (the smooth dotted line in Fig. |Io| ) 
having the appropriate normalisation in this region. 



5.2 Full prior covariance information 

As an indication of the best results one could reasonably ex- 
pect, we now consider the case in which, in addition to the 
azimuthally-averaged power spectra of each input compo- 
nent, we also assume knowledge of the azimuthally-averaged 
cross power spectra between components. These provide 
some limited cross-correlation information between the com- 
ponents at each spherical harmonic mode. Moreover, in this 
case, we also provide the algorithm with the true power spec- 
trum of the CMB input map, with its full acoustic peaks 
structure, rather than the simple functional form assumed 
above. 

Most of the resulting reconstructed maps of components 
are found to be very similar to those obtained in the previous 
section, obtained using no cross-correlation information. We 
therefore do not reproduce them here in full, but list only 
the rms of the residuals maps in Table ^, which we see are 
very similar to those listed in Table ^. 

In this section, we instead concentrate on the few quali- 
tative differences in the reconstructions that occur due to the 
inclusion of cross-correlation information. The most striking 
difference occurs in the reconstruction of the kinetic SZ com- 
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Figure 11. Reconstructed map of the kinetic SZ effect. 
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Figure 12. Power spectrum of reconstructed map of kinetic SZ 
effect (dotted line) compared to the power spectrum of the input 
map (solid line). 



ponent. As shown in the previous section, in the absence of 
cross-correlation information, the kinetic SZ effect is essen- 
tial not recovered at all. Indeed, from Fig. ^, we see that the 
low level emission in the kinetic SZ 'reconstruction' is resid- 
ual contamination from the primordial CMB component. 
If some cross-correlation information is included, however, 
we see in Fig. [ll], that the reconstruction of the kinetic SZ 
somewhat improved. In particular, as found in Paper I, the 
kinetic SZ effect is recovered to reasonable accuracy only in 
those clusters possessing very large thermal SZ effects. This 
clearly shows that the incorporation of cross-correlation in- 
formation between the two SZ effects enhances the kinetic 
SZ reconstruction. The improved quality of the reconstruc- 
tion is also illustrated by its power spectrum, which is shown 
in Fig. |l^ (together with the true power spectrum of the in- 
put map). We see that at least up to I ~ 100, the power in 
the reconstruction matches that of the input map. This is 
a significant improvement on the power spectrum of the re- 
construction obtained using no cross-correlation information 
shown in Fig. B. 

The other major qualitative difference in the reconstruc- 
tion concerns the recovery of the acoustic peaks in the CMB 




Multipole number, / 

Figure 13. The power spectrum of the input CMB map (dashed 
line) compared with the unbiassed estimator (Eu) of the power 
spectrum obtained from the CMB reconstructionfgrey solid line). 
Also plotted are the l-tr error limits. 



power spectrum. In Fig. [13| we plot the unbiassed estimator 
(as discussed above) obtained from the reconstruction. We 
also plot the true power spectrum of the input CMB map 
and the 1-a error limits given by the square-root of the vari- 
ance given in (p7|). To avoid inevitable scatter in the l-cr 
error limits at large I, was approximated by the expo- 
nential function for I >2200. 

Comparing Figs and [y| we observe qualitative dif- 
ferences at high-f . Most importantly, we note that more of 
the acoustic peaks structure present in the true power spec- 
trum is recovered. In particular, the 7th acoustic peak is now 
reproduced with its correct height and position. Beyond the 
the 7th peak, however, error in our simple unbiassed estima- 
tor becomes too large to discern any further detailed struc- 
ture in the power spectrum. This result is in a good agree- 
ment with the analytical predictions by Bouchet & Gispert 
(1999) based on Wiener filtering of Gaussian foregrounds. 

6 DISCUSSION 

In this paper, we have performed an all-sky component sep- 
aration over the whole sky for the Planck mission at full 
angular resolution. The technique employed is an harmonic- 
space maximum-entropy method, which performs a mode- 
by-mode reconstruction of the spherical harmonic coeffi- 
cients of the component maps. Assuming knowledge only 
of the power spectra of the non-CMB components, the algo- 
rithm produces faithful reconstructions of the input maps, 
without the need to apply a Galactic cut. Indeed, it is the 
MEM prior that allows our component separation algorithm 
to cope with such a large dynamics range in pixel values. 
The only component not reconstructed is the kinetic SZ ef- 
fect. We also assumed that the instrumental pixel noise is 
homogeneous and uncorrelated. However, non-homogeneous 
uncorrelated pixel noise can be easily included into the anal- 
ysis. Preliminary investigations of simulations for more real- 
istic Planck observing strategies suggest that the quality of 
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the reconstructions are not significantly reduced, provided 
one obtains reasonable sky coverge. A detailed analysis of 
the effects of scanning strategy on the accuracy of compo- 
nent separation will be presented in a forthcoming paper. 

The power spectra of the reconstructions match those 
of the input maps (expect for the kinetic SZ) on scales 
where the signal dominates the instrumental noise. At high 
£ values, however, the power spectrum of the reconstruc- 
tion consistently underestimates the power spectrum of the 
true map. Indeed, it is easily shown that this estimator is 
biassed. Nevertheless, the component separation algorithm 
also produces an accurate predicted power spectra for the 
reconstruction residuals, which allows the construction of an 
unbiassed power spectrum estimator. This accurately recov- 
ers the input power spectrum at all measured £. 

The effect of including a priori cross-correlation infor- 
mation is also studied. The main effect is to improve some- 
what the recovery of the kinetic SZ effect in clusters pos- 
sessing a large thermal SZ effect. The assumed level of prior 
knowledge of the CMB power spectrum also affects the CMB 
reconstruction, but only at high I. In particular, we find 
that, if ones assumes a CMB power spectrum with no acous- 
tic peaks, the final reconstructed spectrum still recovers the 
oscillatory structure out to the 6th acoustic peak. If, how- 
ever, the true CMB power spectrum is assumed in the com- 
ponent separation process, the recovered spectrum of the 
CMB contains the first 7 acoustic peaks with the correct 
heights and positions. It is also noted that some care must 
be taken in interpreting the recovered CMB power spectrum 
beyond I ~ 2000. In this region, the recovered spectrum is 
dominated by the assumed spectrum. 

The simulations used in this paper did not contain a 
contribution from extragalactic point sources. As shown in 
Hobson et al. (1998), however, the MEM method can accom- 
modate such a components. Moreover, as discussed by Vielva 
et al. (2001), the MEM algorithm can be combined with a 
mexican hat wavelet approach to detecting and removing 
point sources, which leads to a significant reduction in the 
residual point source contamination of the other physical 
components. Recently, Cayon et al. (2001) have presented a 
consistent method for defining the mexican hat wavelet on a 
sphere, which can be combined with the algorithm presented 
here to separate out point source emission in all-sky maps. 
This will be discussed in detail in a forthcoming paper. In a 
future publication, we will also present a detailed statistical 
analysis of the recovery of the kinetic and thermal SZ effects 
in clusters, and discuss its cosmological implications. 
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